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Numerical simulations with access to all possible meson quantum numbers, J PC , are presented 
using two-flavor (up and down) quenched twisted mass lattice QCD with three different lattice 
spacings and four different quark masses. The connection between the quantum numbers (P and C) 
and the symmetries of the twisted mass action is discussed, as is the connection between J and the 
lattice rotation group, for the 400 operators used in this study. Curve fitting of this large data set 
is accomplished by using an evolutionary fitting algorithm. Results are reported for conventional 
and exotic quantum numbers. 

PACS numbers: 12.38.Gc, 14.40.-n, 02.60.Ed 



I. INTRODUCTION 

Although many mesons have been observed in nature, 
the spectrum of known mesons does not include all the 
states expected in Quantum Chromodynamics. Exper- 
imental searches and theoretical studies are continuing, 
but gaps in current knowledge leave room for new ap- 
proaches. Lattice QCD is an established method for ex- 
tracting numerical predictions directly from the underly- 
ing quantum field theory, but lattice explorations of the 
full spectrum of mesons still suffer from limitations. For 
reviews of both theory and experiment for light-quark 
mesons, see Refs. 

A helpful framework for beginning the discussion of 
light-quark mesons is provided by the constituent quark 
model, where mesons are considered to be composed of 
a system of two bound quarks whose spins can couple 
to a singlet (S = 0) or triplet (S = 1) total spin which 
in turn can couple with relative angular momentum L 
between the quarks to produce a total observed angular 
momentum J. Using spectroscopic notation for the states 
2S+1 Lj, one obtains the familiar list of accessible J c , 
as shown in Table Q] 

Since in QCD gluon fields are also present, it is possi- 
ble for gluonic excitations in mesons to contribute non- 
trivially to the observed quantum numbers of a meson. 
Such hybrid mesons with the same J PC as conventional 
mesons are difficult to distinguish but, as shown in Ta- 
ble HI there is a subclass of hybrid mesons, called ex- 
otic mesons, with quantum numbers unattainable in the 
quark model, eg. 0~~, H , 1 h , 2 H , and 3~ + . These 
exotic mesons, for which there is yet no definitive experi- 
mental evidence, would offer a clear signature for excited 
gluon dynamics. 

Twisted mass lattice QCD (tmLQCD) d i] is used 
for the lattice simulations in this work. It is not clear 
a priori whether tmLQCD is a very favorable action for 
lattice simulations of the full meson spectrum. The fact 
that tmLQCD does not respect parity P may be seen as 



TABLE I: States accessible in the constituent quark model 
for different total spin S and orbital angular momentum L 
labeled by spectroscopic notation 2S+1 Lj with corresponding 
J PC . 
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an obstacle, but we will show that it is not insurmount- 
able. Meanwhile, tmLQCD has the advantages of offering 
a cost-effective approach to the chiral limit as well as a 
simple removal of the leading, i.e. O(a), lattice spacing 
errors. In this work, it is shown that tmLQCD's reduced 
symmetries can be understood and accounted for in prac- 
tical simulations, thereby allowing the reader to make an 
informed decision as to whether this option is preferable 
for future studies of the meson spectrum. The issues spe- 
cific to tmLQCD that affect our operators are discussed 
111 Section HTT1 

Much of the discussion of meson operators is not 
specific to tmLQCD. See Section [XT] for this tmLQCD- 
independent presentation. In particular, we choose me- 
son correlators that can be built from quark and anti- 
quark propagators originating at a single lattice site aug- 
mented by gauge field links defined on extended spatial 
paths. 

All possible quantum numbers are obtained by using a 
variety of options for the paths of gauge fields connect- 
ing quark to anti-quark. (Gauge fields are subsequently 
smeared at both source and sink, but quarks are only 
smeared at the sink in this work.) With just one quark 
propagator inversion, these "excited glue" operators are 
minimally expensive, and the results of our simulations 
allow us to tabulate the relative strengths of the overlaps 
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that various mesons have with various operators. It is in- 
teresting to see the extent to which these "excited glue" 
operators couple to exotic mesons, and also to conven- 
tional quantum numbers. 

Meson masses are extracted from lattice QCD simula- 
tions by fitting a linear combination of exponentials to 
correlators as a function of Euclidean time. It is a deli- 
cate business. For example, the number of exponentials 
that should be used for a certain fit depends on the par- 
ticular channel being studied and also on the quality and 
quantity of data. Our study of the meson spectrum will 
include hundreds of correlators, some of which should be 
fit simultaneously since each meson will appear in mul- 
tiple correlators. Because sink smearing is more easily 
varied than source smearing due to the expense of re- 
computing quark propagators, we want a fitting method 
that does not require computation of a complete corre- 
lator matrix, i.e. we want the freedom to consider more 
sink options than source options. 

To address all of the delicate issues of fitting, the evolu- 
tionary algorithm introduced in Refs. @, 0] is used. The 
workings of the algorithm are well-understood in terms of 
the basic principles underlying biological evolution, but 
it is a black box algorithm in the sense that human in- 
tervention, and therefore human bias, is avoided. For 
example, the algorithm will identify the number of expo- 
nentials that minimize the x 2 / n dof for a given fit. This 
black box method is also general enough to handle multi- 
correlator fits with no need of a complete correlator ma- 
trix. The data-fitting technique proposed in Refs. @, @] 
is independent of tmLQCD. This is its first application 
to such a large set of lattice QCD data, and the imple- 
mentation is described in Section [IVI 

Section [V] explains the parameter choices used in our 
numerical simulations. The present study is exploratory, 
and we have therefore chosen to perform quenched sim- 
ulations. 1 Lessons learned in this study about extended 
meson operators, fitting algorithms, tmLQCD and the 
meson spectrum are applicable to future studies beyond 
the quenched approximation. Section IVII presents and 
discusses the results. 

Section IVIII draws some conclusions about the meson 
spectrum. It also highlights properties of our chosen op- 
erators, comments on the appropriateness of the twisted 
mass action to this physics, and underscores the valuable 
qualities of the evolutionary fitting technique. 



II. LATTICE-SYMMETRIZED MESON 
OPERATORS 

Lighter quark mass calculations require improved 
statistics for hadron mass resolution. This requirement 



1 For the status of dynamical tmLQCD, see Ref. Q and references 
therein. 



TABLE II: Number of copies of the irrep A of O in the 
reduction of the subduced representation of the continuum 
irrep J of SO(3). 
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and the need to disentangle physical states with the same 
quantum numbers may be accomplished in part through 
the creation of more operators that represent the chan- 
nel in question. The creation of more elaborate operators 
also allows for the study of hybrid and exotic mesons. 

Operators with displaced quarks have been con- 
structed using group theoretical techniques but their us- 
age requires the calculation of quark propagators from 
multiple lattice sites. (See, for example, Refs. [9l. [Tol. [TT| .) 
Consideration of the operators available through gluonic 
extension alone therefore is numerically expedient. More- 
over, it allows investigation of the coupling of operators 
with "excited glue" to conventional and exotic states. 
This section provides a complete discussion of the opera- 
tors (first introduced in Ref. [12]) used in our simulation. 



A. Lattice symmetry group 

While parity (P) and charge conjugation (C) may be 
conserved by lattice actions, the continuous rotational 
symmetry of nature is broken and one requires opera- 
tors adapted to the symmetry group of the lattice. For 
mesons this is the octahedral group O with 24 elements. 
The group O has five conjugacy classes conventionally la- 
beled {E, 3C| , 8C3, 6C4, 6C2} and therefore admits five 
(unitary) irreducible representations (irreps): two one- 
dimensional irreps Ax and A 2 , one two-dimensional irrep 
E, and two three-dimensional irreps T\ and T 2 . The di- 
rect product of the parity, charge conjugation, and octa- 
hedral groups is denoted O pc . 

For an operator adapted to the representation A PC , 
where A G {Ai, A 2 , E, Tx, T 2 } is an irrep of O and 
P,C € { + ,—}, one identifies the possible physical states 
J PC to which it corresponds using Table |TT] which 
shows the number of copies n]( of irrep A contained 
in the reduction of the subduced continuum rotation 
group (SO(3)) irrep J [13j . For example, an operator 
transforming as the irrep E~* could have spin content 
(J PC — 2^ , 4 H , 5 H , . . .). In theory one needs many 
operators of each A PC to resolve the range of physical 
spins J PC in the tower of states to which A PC corre- 
sponds. 

Table [TT| is arrived at through a consideration of the 
character table of O shown in Table [TTT1 The trace of the 
rotation group matrix for spin J and rotation angle 
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TABLE III: Character table of O showing the character \ 
for the given class £ in the irrep A. The final line displays 
the angle of rotation 9 common to the elements in each class 
required in Equation (TTJ. 
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FIG. 1: Diagrammatic representation of the building block 
operator. 



corresponding to class £ is given by 2 

X {J) (0 = sin[(J + 1/2)0(0]/ sin[0(O/2] . (1) 

Those matrices in the continuous irrep J which corre- 
spond to the subgroup O of 50(3) form a representation 
of O which is, in general, now reducible. Equation ([1]) 
combined with Table IIII) allows one to determine the 
subduction of continuum J to discrete A of O via the 
decomposition formula 

ni = -Y / PiX iA) (0*X {J) (0 (2) 
9o ^ 

for the number nj^ of copies of A in the subduction of J. 
Here p^ is the number of elements in class £ and g is the 
order of the octahedral group. 



B. Operator building blocks 

One may construct zero-momentum operators trans- 
forming as irreps of O from the space of operators 
spanned by 

M jM {t) = y £^ a {x)U j , k {x)<4i h {x) , (3) 

X 

where the gauge link part of our operators is defined via 

U jik (x) = U^Ukix+^U-jix+j + kp-kix + k) . (4) 

Here j, k = ±l,±2,±3;j 7^ k and a and b are spinor 
indices for a total of 24 x 16 = 384 operators. The j de- 
notes a four- vector of unit length along the spatial axis j. 
See Figure [1] for a diagrammatic representation of the 
building block operator. Superpositions of such opera- 
tors for meson spectrum analysis have been sug gested 
and used previously in special cases 

BUI E! Eg. Eq- 
uation of correlators constructed from operators of this 



2 This formula is to be interpreted as a limit in the event the 
denominator vanishes. 



form requires only calculation of propagators from a sin- 
gle source. Construction of operators transforming as 
irrep A PC is facilitated by observing that one can effec- 
tively consider the transformation of the spinor and link 
paths independently and then combine them via octahe- 
dral Clebsch-Gordan coefficients. Link smearing can be 
performed on operators at both the source and sink of 
the correlator, and quark smearing at one end to sup- 
press high energy states with no change to the symmetry 
of the symmetrized operators. 

C. Spinor contribution to group structure 

The contribution to the group structure due to 
spinor indices is determined by the 16 bilinears ipFip, 
where F represents one of sixteen 4x4 matrices, 
{/, 75 , 74 , 7475 , li , lilb , 04* , £y 7c Cjfc } • The first four bilin- 
ears are scalars while the last four three-index objects 
are vectors under the rotation group in the Euclidean 
continuum. 

Parity (P) and charge conjugation (C) of the bilinears 
are identified via 

C^=^f, Vfl>i=wl>, (5) 

c^t = - (cy>) T , = V74, 

where C is the matrix implementing charge conjugation 
and we have suppressed the action on coordinates, which 
under parity sees x = {x\, X2, £3, X4) — > Xp, where 

x P = (-x 1} -x 2 , -x 3 ,x 4 ). (6) 

To classify the bilinears into irreps of O one uses its 
character table found in Tabic iLLTl to project out the oc- 
tahedral irreps from the representations generated by the 
bilinears. Since the bilinears are scalars and vectors in 
the continuum it follows that their respective spans are 
also invariant under O. By inspection it is straightfor- 
ward to show that the character table for the scalar and 
vector bilinears is identical with that of A\ and T\ re- 
spectively. Since the multiplicity n J A of irrep A in any 
representation J is given by Equation @, it follows triv- 
ially that the scalar and vector bilinears form the basis 
of A\ and T± irreps respectively. 

Furthermore it may be verified that the i th vector bilin- 
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TABLE IV: Matrices F for the generators Ci y and d z of the 
octahedral group for each irrep A. 
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FIG. 2: PC-symmetrized basis elements for gauge field links. 
We have suppressed the overall normalization factor of 1/2. 



TABLE VI: Character table showing x (U '(0 for represen- 
tation U PC induced by UfF . Here positive and negative are 
to be interpreted as +1 and —1 in the entries. 

E 3C'i 8C 3 6C 4 



6C 2 



6 2P (P + 1)C 



D. Link contribution to group structure 



TABLE V: Octahedral symmetries of the spinor bilinears. 
The operators have been separated into quadrants according 
to their dependency on twist angle as discussed in Ref. [l2| . 
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ear component transforms as the i row for our choice of 
matrix representation r^ Tl ^(i?) given in Table [TVl Since 
the action on the spinor under rotation R is ip — > ipS(R) 
and ip — > this amounts to verifying that 

S(R)F i S^(R)=J2^ l) (m, (7) 

3 

for the elements R of O and the four different vector bi- 
linears F. Here S(R) is the matrix implementing the 
rotation on the spinors. It is sufficient to verify Equa- 
tion (O for the two generators, namely the 7r/2 rotations 
about the y and z axes, C± v and C± z , for which these 
matrices are given by 

S(C± y ) - -±=(1 + 7173), (8) 

S(C iz ) = -±=(l + 7 27i)- (9) 

Hence the reduction of the spinor structure of our op- 
erators is given in Table |Vl where now each bilinear may 
be classified uniquely by its irrep A PC , row A, and irrep 
multiplicity index a as F^ 



As with the spinors, one can simplify the discussion of 
the rotational properties of the gauge links by consider- 
ing first the parity and charge conjugation and only then 
the rotations in O. Parity of a link about the point x 
is implemented via inversion about all three spatial axes 
as usual. Charge conjugation sees U — * U*, however 
to compensate for an overall transpose arising from the 
transformation of our quark bilinear, we formally trans- 
form U — > U* to achieve the correct overall charge conju- 
gation properties of our operators. The parity and charge 
conjugation contributions due to this link structure can 
then be taken into account by defining the PC-adapted 
superpositions, 

u if = \ ( U i* + p U- 3 .-k + CU ktj + PCU-^-i) , 

( 10 ) 

where k is defined via k = i x j. Diagrammatically, 
Vff is shown in Figure O For fixed P and C the space 
spanned by Ufj is invariant under O and will gener- 
ate a representation U PC of the group. A basis for 
the six-dimensional space may be found by restricting 
(*,i) to {(1, 2), (1,3), (2, 3), (2,1), (3,1), (3, 2)} and this 
will in turn induce the particular matrix representation 

r (c/PC) (R) 

In order to determine how each of these four six- 
dimensional representations reduces, it is convenient to 
calculate their respective character tables. A straightfor- 
ward consideration of the effect of a single element from 
each class of O on the basis vectors allows one to find 
the character table for U PC shown in Table IVI1 Use of 
Equation ^ with J as U PC reduces the latter to irreps 
A PC of O pc as shown in Table IVIII So, for example, 
U++ = A++ © E++ ® T++. 

To arrive at the explicit form of the lattice- 
symmetrized gauge links, note that, by Table IVIII each 
of our spaces reduces to no more than a single copy of 
each irrep. For such simply reducible representations the 



TABLE VII: Reduction of U PC to irreps A PC of O 1 
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reduction is straightforward 3 and in our case one has the 
following formula for the symmetrized link fields: 



(x) 



(11) 



ReO ^(i,j)(m 



SO. T> 



(m,n)(m,n) W 1 MM W 



where (m, n) and are chosen so the denominator does 
not vanish and the accessible A PC are taken from Ta- 
ble IVII1 Here T^J (R) is the matrix representation for 
the irrep A of dimension <Ia taken from Table IIVI Low- 
ercase Greek denotes row indices of the group O, and 
the sum is taken over the six pairs listed above. 

The {/■ (x) transform as the row A of irrep A PC and 
may be uniquely identified by these latter parameters. 
(Since each A PC occurs only once in Table IVIII there is 
no need, unlike in the spinor case, to further identify a 
multiplicity for the irrep.) Using the explicit irrep gen- 
erator matrices in Table IIVI to construct the full matrix 
representations for each irrep, one can evaluate Equa- 
tion (fTTj) to produce the final lattice-symmetrized gauge 
links. The results are tabulated in Table [VTTTl A diagram 
of a lattice-symmetrized gauge link structure is shown in 
Figure [3] 

In Equation (fT0|) , all the link structures with PC ^ ++ 
are suppressed by powers of lattice spacing relative to 
PC = ++. Further factors of lattice spacing can appear 
when Equation (jTUJ) is used to form a gauge structure like 
the one shown in Figure [3l Such overall factors of lattice 
spacing will not be of direct relevance to our study of 
the mass spectrum, but would be of greater interest for 
decay constants and matrix elements. 



E. Total representation reduction and operator 
construction 



Having classified the spinor and link components of 
our operators into irreps of O p , it remains to combine 
them and reduce all possible direct product representa- 



3 See Ref. [lTl page 74] where the Clebsch-Gordan coefficients for a 
simply reducible direct product representation are given. In the 
case at hand one simply considers an arbitrary representation in 
place of the direct product to arrive at our result. 
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FIG. 3: Extended gauge structure U x 
Parity and charge conjugation symmetries are readily verified. 
The Ti irrep has three rows, this being the A = 3 component 
which is symmetric about the z-axis. Illustrations of other 
lattice-symmetrized gauge structures may be found in Ref. 0]- 



tions. Parity and charge conjugation for the product rep- 
resentations are given by 



P = P f P u c = c f c u 



(12) 



where / corresponds to the fermionic part and u cor- 
responds to the gauge part. Noting that the character 
of a direct product representation like A / ® A u satisfies 
X (A / ®A„)(^) = x {Af) {hx iAu) (0 and usin S Equation © 
with J = Af (g) A u and the characters in Table IIII| one 
may reduce the octahedral part of each direct product 
into irreps of O as shown in Table IIXI 



The results shown in Tables fVl IVIII and IIXI along with 
Equation (fL2|l finally allow the 384-dimensional represen- 
tation generated by the space spanned by Mj^, a ,b{t) to 
be reduced into irreps of O pc as shown in Table [Xj It 
is of note that while the local quark operators can ac- 
cess only irreps of the form A PC and T PC of PC , the 
addition of extended gauge field structure admits opera- 
tors of every possible irrep. In principle, then, this set of 
operators can couple to all meson states. 



The construction of the lattice-symmetrized operators 
themselves is accomplished by combining the spinor and 
link operator components F A and using the 

same formulae for parity and charge conjugation given 
in Equation (| 1 2j) and by combining the irreps of O using 
Clebsch-Gordan (C-G) coefficients for the reduction of a 
direct product of two irreps of O, Ai g) A2. As suggested 
by Table HXl the direct product representations for irreps 
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TABLE VIII: Shown are the lattice-symmetrized gauge links, U\ (x), derived from Equation in terms of the PC-adapted 
gauge field structures Ufj(x). Each subtable corresponds to a different PC combination and contains the lattice-symmetrized 
gauge link combinations along the top, denoted by A A , and the PC-adapted structures along the left, denoted by U(i,j). For 
irreps of dimension one the superscript denoting the row A of the lattice-symmetrized link structure is suppressed. The actual 
coefficient equals the sign of the entry times the square root of the absolute value of the entry in the table. Thus, for example, 

one has [7 3 1 = (yi~2 + + ^23 — ^21) /2, which in turn can be expanded via Equation (|10fl to arrive at Figure [3] 
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TABLE IX: Reduction of product representations A/ ® A„ 
into irreps of O. Only irreps A\ and Ti are accessible for 
Af. Reduction of other octahedral direct products may be 
inferred from Table KU 
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have been calculated and are given in Table IXII The 
lattice-symmetrized operators, (t), are finally 
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TABLE X: Reduction of the representation generated by the 
span of Mj,k,a,b(t) into PC irreps. The 384-dimensional rep- 
resentation reduces into the sum of the irreps listed in this 
table with multiplicities shown. 
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of O are simply reducible, allowing one to use the formula 
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to determine the C-G coefficients. (See 
Ref. [13, page 74].) Here //, fii, and [i 2 are chosen 
so the denominator does not vanish. Using the irrep 
generators given in Table IIV1 these C-G coefficients 



where the allowed irreps for the spinor bilinear and link 
components are determined by Tables |V] and IVIII These 
fix P and C for the operator while the irreps A of O 
are those allowed by the C-G series in Table IIXI Dirac 
indices on F and color indices on U and both indices 
on the spinors have been suppressed. We note that each 
operator is thus uniquely identified by its irrep A, its 
row A, and the direct product from which it originates, 



A 



f 



III. TWISTED MASS 

It is a feature of twisted mass lattice QCD (tmLQCD) 
that, at maximal twist, 0{a) errors are absent from phys- 
ical observables 0, IE). Critical slow-down is softened 
in tmLQCD and thus it permits simulations with light 
quark masses. Unlike theories that require the tuning of 
separate parameters for the improvement of each individ- 
ual operator, tmLQCD requires the tuning of only one, 
the standard mass parameter tuq. At maximal twist, in- 
formation about the physical quark mass is given by the 
twisted mass parameter /iq. However, the theory modi- 
fies the parity symmetry of QCD. The implication of this 
for meson correlators is examined in this section. 
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TABLE XI: Octahedral group Clebsch-Gordan coefficients derived from Equation (|13p . The actual coefficient equals the sign 
of the entry times the square root of its absolute value. Superscripts on irreps denote the row. If the irrep is one-dimensional 
the superscript is suppressed. Since the C-G coefficients for Ai ® A2 are the same as for A2 (g> Ai only one combination is listed. 
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The tmLQCD fermion action is 

S F = a i ^2Q(x)S~ 1 (x,y)q(y) , 

q=u,d x,y 



where 

(15) ^2S- 1 {x,y)q(y) 

V 

= ^2 7p [Un(x)q{x + a/t) - Ufa - afi)q(x - aft)] 
^~ [ u v( x )l( x + a A) + Ufa - ajj,)q(x - a/t) - 2q(x)] 



[m + ijsVq] <l{x) 



(16) 
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and 



(17) 



is the twisted mass parameter. The tmLQCD action 
of Equations (|15H16p is displayed in the "twisted basis" . 
Fermion propagators are computed in that basis and then 
converted to the "physical basis" via 



(»5Vi)physical 
(Geophysical 



(S M )tw is te d e^ /4 , (18) 



_ „-i-y5w/4 



{Sd) twisted 6 



Except for Equations (|15H16p and (|34|) . all discussions in 
this work are in the physical basis at maximal twist. 

While parity P defined in standard fashion by Equa- 
tion |5]) is not a symmetry of the tmLQCD action, flavor- 
parity 



P = FoP 



(20) 



which combines the parity operation with a 180° rotation 
about the second isospin axis, 



(21) 



= — d, Ti<IT\ = u, 
J-2u!F\ = — d, T^dT\ = u, 

is respected In our notation, JF 2 is defined by 

J~2 = e l7rra (22) 

where r 2 is an isospin operator, defined in Ref. [ic[ . whose 
action on quark fields can be inferred from the commuta- 
tion relations in Equations (31) and (32) found therein. 

Notice the relations involving reflections, 

S- 1 (0,y ) U) = ia 5Sf(0,yr,U I ) l5li , (23) 

where the subscript / refers to the lattice coordinates and 
links after inversion in the ith direction. Applying this 
relation in all three spatial directions and a consideration 
of the inverse leads to a flavor-parity relation 



S u (x,y;U) = jiS d (xp,yp; U P )-fi 



(24) 



where xp is defined by Equation (J6j) and Up represents 
the links after inversion under spatial parity. The tm- 
LQCD action also preserves charge conjugation, 

S U (U) = CSl{U*)& (25) 

and has a Hermiticity relation, 

S U (U) =75^(C/) 75 . (26) 

In our simulations we consider charged-meson two- 



point correlators, 
Cab® = 



(27) 



^2 u(0) l4 (F A U A (0)) t - /4 d(0)d{x)F B U B (x)u(x) 



\ x, spins 
colors 



and neutral-meson two-point correlators with discon- 
nected contractions omitted, 



(19) N AB (t) = 



u(0)-f 4 (F A U A (0))^ 4 u(0)u(x)F B U B (x)u(x) 



\x, spins 
colors 

~(U <r-> (J) , 



(28) 



where each of U A and U B is one of the symmetrized link 
fields defined in Equation (fTTjl . and each of F 1 

and F J is one of the 16 unique products of Dirac matri- 
ces FA defined in Table fVl Choosing F A to be Her- 
mitian or anti-Hermitian and evaluating the (connected) 
contractions, a few lines of algebra leads to valuable con- 
clusions: 

• In the configuration average, Caa® is real for 
tmLQCD fermions. 

• In the configuration average, Naa® is real for 
tmLQCD fermions. 

Note that the imaginary part of Naa® does not vanish 
if the flavor symmetrization, (u <-> d), is omitted. In 
practice, correlators in the charged case are averaged as 
well but for improved statistics. 

We are interested in correlators Cjj(t) and Njj(t) be- 
tween the fully-symmetrized source and sink meson oper- 
ators / and J of the form given by Equation (|14|l . 
These correlators involve C-G superpositions of Cab® 
and Nab® respectively and the above results hold for 
C'jjit) and N' IT (t) as well. 

It is necessary to consider what effect the use of a 
twisted mass action will have on the group theory dis- 
cussion from Section[Tl] The distinction between respect- 
ing P and not P is clarified by using the three symme- 
try relations, Equations (|24M26|) . to evaluate two-point 
correlators of operators of given parity P. Table IXIII 
shows the resulting orthogonalities among various cre- 
ation/annihilation operators. Note that charged mesons, 
containing ud or du creation operators, have different re- 
lations from neutral mesons, containing uu or dd creation 
operators. 

In Table IXII| the charged meson entries that vanish 
without a flavor interchange (u <-» d) are direct conse- 
quences of PC symmetry. This is equivalent to the phys- 
ical PG symmetry since G = F^C. These lines in the 
table prove that PG = — 1 states (i.e. A ++ and A ) 
are orthogonal to PG = +1 states (i.e. A" 1 and A '"). 
For example, because the physical pion is an eigenstate 
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TABLE XII: Orthogonality relations among creation/annihilation operators in tmLQCD. This table applies to any fixed 
A 6 {Ai, A2, E, Ti, T2}. Superscripts denote P (not P) and C for the operators / and J. The charged mesons are only 
eigenstates of G-parity with eigenvalue G = C( — l) 1 = —C, but we abuse the notation here and elsewhere in the paper. 



A^ 

A++ 
A++ 
A++ 



J 



charged meson relations 



{^iP^iSuJSd) - (u <-> d) = 
(jilhAJSd) = 
(74l t 74S„ JS d ) = 

(74^74^ + (u <-> d) = 



neutral meson relations 
{-fil^liSuJSu) - (u <-y d) = 
{'YA^liSuJSu) = 
(^d^liSuJSu) + (u« d) = 
^/il^iSuJSu) = 



A+- 
A+- 
A+- 



A+^ 
A+- 
A^ 
A— 



{-fil^iSuJSd) = 
jil^^SuJSd) - (u d) = 
j4l if l4S ll JSd) + (u <-> d) = 



(^d^liSuJSu) - (u <h> d) = 



A+^ 
A+- 
A 1 



A+- 


A— 


(7 4 / t 745'u^5 , d) 


= 


("fil^liSuJSu) - 


■f (it <-> d) 


= 


A"+ 


A++ 




= 




■f (u «-> d) 


= 



A-+ 
A-+ 
A-+ 



A" 1 
A- 
A" 



7 4 / t 74<5 , M^<S'd) + (" <-> d) = 
j4,Pj4S ll JSd) - (u f-» d) = 
(74^74^ JSrf) = 



(j^jiSuJSd) + {u <-> d) = 
(liPliS u JS d ) = 
{liPliS u JS d ) = 

(jil^iSuJ Sd) - (u <-> d) = 



(74/ t 745 , M dS'n) = 
{lil^liSuJSu) — (u d) = 
("f^jiSuJ 'S u ) = 



<74/ T 74^u J5 U ) = 
{lil^liSuJSu) + (u <-» d) = 

(74^74^ J5 U > - 
{lil^liSuJSu) — (it <-> d) = 



A" 
A" 
A" 
A" 



A+^ 
A+- 

A— 



of PG with eigenvalue +1, we know that no charged pion 
signal can appear within Af + nor within A[ . 

Unfortunately, the charged physical pion will in gen- 
eral couple to both the A\ and A^ + channels. This 
is because it is an eigenstate of P and G separately, 
whereas tmLQCD does not respect these as good quan- 
tum numbers. Instead, tmLQCD respects P and C, but 
the charged pion is not an eigenstate of those. Notice that 
the lines in Table IXLTl that involve «-> d)" merely 

show that unphysical states of opposite flavor-parity are 
orthogonal, i.e. 

(^il^iSuJSd) + (u <-> d) 
= (drfiP^uuJd') + (uj^PjiddJii) 
= ( (d-f 4 P 74U + -u7 4 / t 74 d) (uJd + dJu)) . (29) 

Each factor in parentheses has definite flavor-parity, so 
finding that this matrix element vanishes (as was found in 
some entries in the table) means that the flavor-parities 
are orthogonal. Similarly, the charged physical mesons 
in and A\ cannot be separated either. 

In Table IXII1 the neutral meson entries that vanish 
without using (u <-> d) are direct consequences of C sym- 
metry. These lines prove that C = +1 states (i.e. A ++ 
and A h ) are orthogonal to C = — 1 states (i.e. A H and 
A ). For example, because the physical neutral pion 
is an eigenstate of C with eigenvalue +1, we know that 
no neutral pion signal can appear within A\ nor within 

In contrast to the case of the charged pion, the phys- 
ical neutral pion is an eigenstate of P and since tm- 
LQCD respects P the pion will couple to A\ + and not to 



A'l . However, the AJ + channel will couple to the flavor- 
singlet pseudoscalar (i.e. the SU(2) if) and if we neglect 
disconnected diagrams, then the mass of the SU (2) tj' is 
identical to the mass of the neutral pion. This degener- 
acy is not specific to the pion-?/ system; it occurs for any 
angular momentum irrep with any C . Therefore A PC 
contains the same spectrum of masses for both P values 
in neutral channels. 

As discussed in Ref. [12] , all A PC combinations can be 
obtained from operators that are independent of twist 
angle. Since the present work is restricted to maximal 
twist, those operators are not emphasized. In cither 
case, tmLQCD simulations will still contain mixing as 
discussed above: For charged mesons we cannot separate 
the A ++ and A~~ pair, nor can we separate the A^ 
and A ^ pair. For neutral mesons we cannot separate 
the A ++ and A ^ pair, nor can we separate the A^ and 
A pair. A striking example is the appearance of a pion 
signal for the operator u(x)~fid(x). 



IV. CURVE FITTING WITH EVOLUTIONARY 
ALGORITHMS 

There are several motivations for our use of evolution- 
ary fitting. In general, the sheer quantity of fitting to be 
done requires an automated black box method, an objec- 
tive being pursued by others as well [18| . An important 
part of this goal is that the method not depend on any 
subjective parameters such as timestep fit ranges, a fixed 
number of states, or initial parameter values. This is 
not simply to speed up the fitting process but for re- 
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producibility as well. While evolutionary fitting requires 
mutation and breeding steps to be identified as well as 
probabilities to be fixed in the algorithm, the goal of 
minimizing the x 2 / n dof is objective. Any two fit algo- 
rithms, evolutionary or not, can have their results on the 
same data readily compared by this statistic and either 
have the fits confirmed to be equivalent or the better fit 
selected. 4 

Other reasons for using this fitting method are specific 
to our problem. We require a method that extracts a set 
of common states from multiple correlators in the same 
channel. This is often accomplished using the variational 
fitting method [H, HU HH- Due to the large number 
of operators present, however, we run only the diago- 
nal correlators of the correlator matrix and not the full 
correlator matrix generated by putting all possible com- 
binations of operators in a given channel at source and 
sink. The latter would have been required to use the 
variational method. We also smear our source and sink 
operators differently, which results in a correlator matrix 
which is not Hermitian, which is also required for that 
method. Evolutionary fitting in principle allows us to 
identify continuum angular momentum states (J) across 
lattice symmetry channels which have only octahedral 
symmetry. Finally, a method that is able to identify po- 
tentially weak contamination signals occurring due to the 
use of twisted mass is also of value. 

We have detailed the algorithm itself in Refs. @ 
and Q. Here a summary is presented to provide con- 
text for evaluation of the method. 

The idea behind an evolutionary algorithm is to con- 
sider candidate solutions to a problem as individual or- 
ganisms in a population. This population is allowed to 
mutate and breed to produce successive generations. The 
types of operations involved in mutating an individual 
or breeding a pair of individuals will be specific to the 
problem at hand. These must be sufficient for the organ- 
isms to have a good probability of being able to explore 
the solution space. Coupled with a function to measure 
the fitness of an organism which in turn influences which 
organisms will be allowed to populate the next genera- 
tion, evolution is allowed to take its course with succes- 
sive generations approaching a better solution. Different 
individuals in the population will have desirable charac- 
teristics which will be disseminated over time with high 
probability to the rest of the population. In a sense, the 
population is able to explore the solution space to the 
problem in parallel, making it an effective technique in 
general for finding good, and sometimes unexpected, so- 
lutions to complex problems. 

In the context of lattice QCD, the organism is a fit 
function that is a linear superposition of an a priori un- 
known number of exponential functions whose exponents 



4 With such comparisons, fitting algorithms themselves may 
thereby evolve. 



determine the spectrum of masses in the correlator. The 
simplest case involves a fit to a single correlator and we 
have detailed this elsewhere [1, 0|- In the current situ- 
ation we have multiple correlators for a single channel 
sharing the same energy states {E m : m = 1, . . . , m max }. 
The total fit function to all the data is then a set of fit 
functions, one per correlator, each of the form: 

G«(i) = J2 Z( n l) + e"^ i)(T_t) ) . (30) 

71=0 

Here i is the correlator index, Umax is the (variable) num- 
ber of states found in that correlator, Z n is the coefficient 
for the energy state E] n , and T is the temporal extent of 
the lattice. The form of the fit function reflects the fact 
that we are only fitting diagonal correlators in the cor- 
relator matrix (essentially the same operator at source 
and sink though we do allow for different smearings at 
source and sink), and that we are using periodic bound- 
ary conditions. The desired best fit minimizes x 2 / n dof 
of the total fit function G. Assuming the datasets are 
not correlated in any way, % 2 (G) is simply the sum of 
X 2 (G^), the correlated x 2 on t ne i th correlator [7]. The 
number of degrees of freedom is 

ndof{G) = n data - m max - ^ n^ ax , (31) 

i=l 

where ndata is the number of timesteps fit in each corre- 
lator times the number of correlators (i m ax)- The fitness 
of the organism, /(G), which we desire to maximize is 
therefore defined to be — x 2 (G)/ndo/(G). One notes that 
the degrees of freedom fluctuate with the number of pa- 
rameters (states and coefficients) in a particular organ- 
ism. The complication of searching such a discontinuous 
function space as well as the independence from initial 
conditions (the initial population is chosen at random) 
are some of the principal advantages of the evolutionary 
fitting method. 



The information required to construct a given fit or- 
ganism is coded in its genotype. The subfit for dataset 
number i is represented by a list of Umax coefficients 
(Zn\ In ) where I € {1, . . . , m ma x} is an integer index 
indicating to which of the m ma x energy states E m the 
coefficient is associated. In summary, for a fit of i ma x 
correlators the complete genotype is of the form: 

Fit Genotype 
= (Dataset coefficients, Mass list) 
= ((Dataset 1 coeffs, . . . , Dataset i max coeffs), 
Mass list) (32) 
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with 

Dataset i coeffs 
Mass list 



• (Z% ,1% )) 

VI ' ' VI K ' 



(33) 



With the above genotype defined, operations for mu- 
tation and breeding become apparent. The genotype is a 
hierarchy of lists so we coded procedures that proceeded 
recursively through the structures present. 5 Mutations 
of lists involve mutating the elements in the list. For a 
list of elements of the same type, adding or removing a 
random element, and if order is meaningful to the list, 
reordering it, are other possible mutations. Breeding (or 
crossover) of two lists will produce two new lists con- 
taining parts of each, and potentially of different lengths 
if the lists contain elements all of one type. Ultimately 
one has to mutate numbers, either reals or integers, and 
this can be accomplished by adding a Gaussian random 
variable or flipping bits in a binary representation respec- 
tively. (Our integer index has to be interpreted modulo 
the number of energy states m max so that it always maps 
to an individual state.) Breeding of numbers can be done 
by randomly interpolating the real numbers or exchang- 
ing subsets of the bits of the integers. 

It is of value to introduce special mutations to the full 
genotype. One mutation does a fixed number of steps of a 
Newtonian optimization of the fit function of a single or- 
ganism. Here we used the Levenberg-Marquardt [23, HH 
method with the current functional form defined by the 
organism's genotype and the values of its individual pa- 
rameters as the initial conditions. Also we introduced a 
reduction mutation which would convert genotypes that 
represented the same function into a common form so 
that the fitting algorithm would converge to a single rep- 
resentation of the solution. Consult Refs. [6|,|7| for further 
detail of these steps. 



V. SIMULATION DETAILS 

The quenched configurations used for this work were 
previously discussed in Ref. 25], and the details are re- 
produced in Table IXIIII of this work for completeness. 
In addition to tmLQCD fermions, a Wilson fermion was 
considered for purposes of comparison and those parame- 
ters are also listed in Table IXTITl Note that at each of the 
three (3 values, the four quark masses are comparable by 
their approximate ratios with the physical strange quark 
mass: m s , m s /2, m s /3, and m s /6. The Wilson quark at 
(3 = 6.0 was chosen close to m s /2. 

The tuning to maximal twist was performed by varying 
mo for each p,Q until cu became ir/2 as defined by [26l. 1271. 



tancj 



(34) 



where P, U M , are the local bilinears for charged mesons 
with pseudoscalar, vector and axial vector quantum num- 
bers respectively. 

The diagonal correlators of the 384 extended and 16 
local operators from Section HT1 were calculated for each 
lattice spacing using degenerate quarks at our four quark 
masses in both the charged and neutral channels for 
twisted mass and at the single quark mass for Wilson. 
Neutral twisted mass channels did not include the so- 
called disconnected contributions, i.e. contractions be- 
tween a quark and anti-quark within the source (or sink) 
operator. 

Smearing was used to reduce contributions from ex- 
cited states to correlators at the shortest time extension. 
Gaussian quark smearing was performed at the sink and 
stout link smearing at both source and sink as described 
in Ref. [2{|. Quark smearing parameters were a = 0.15 
and n a — 64 for all lattice spacings. Stout link smear- 
ing used p = 0.15 for f3 = 5.85 and p = 0.2 for /3 = 6.0 
and 6.2 and n p = 16 for all three lattice spacings. 6 For 
comparison, correlators of unsmeared operators were also 
computed. 

Before fitting, diagonal correlators corresponding to 
operators that differed only in their row A were aver- 
aged since these must be the same statistically by symme- 
try. Specifically, the diagonal correlators from operators 
in Equation flH) with the same (A PC , A" fPfCf , A^"* 7 ") 
were averaged, leaving the number of correlators to be fit 
for each channel given by Table IXl for a fixed quark mass, 
fixed type (twisted mass charged, twisted mass neutral, 
and Wilson), fixed lattice spacing, and fixed smearing. 

For each such set of correlators, evolutionary fits were 
done to the data at all timesteps using an overall pop- 
ulation of 480 organisms. 7 Two fits were done for each 
dataset to test fit consistency. The first was stopped at 
exactly 600 generations. The second run, the results of 
which were used, performed at least 600 generations but 
was allowed to continue up to 1200 generations, stopping 
in between only if no improvement in the best genotype 
of a given generation was seen for 200 generations. The 
genotype was allowed to contain up to 8 masses, and 
the fit for each correlator could have up to 8 coefficients 
pointing to elements of the mass list. All coefficients and 
masses were restricted to be positive. Once the evolu- 
tionary algorithm had obtained a best fit function for the 
dataset, bootstrap errors [3(| were generated by fitting 
this function to bootstrap configurations. 8 



5 A modern list-based language, Python [22l |. was convenient for 
this implementation and fast enough for our purposes. 



6 See Ref. _29] for explicit definitions of these parameters. 

7 The population was distributed over 4 islands with parameters 
Neiite = 5, N diversity = 5, Nmutant = 20. See Refs. [EI3- 

8 Bootstrap fits were done using the Levenberg-Marquardt 
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TABLE XIII: The parameters used for simulations in this work. Lattice spacings are taken from Ref. [26| using ro = 0.5 fm. 
Each (amo,a/io) pair is the result of tuning to maximal twist as discussed in Ref. [2y|, except the Wilson case of course. 





a [fm] 


#sites 


^configurations 


amo 




twist angle (degrees) 


5.85 


0.123 


20 a x 40 


600 


-0.8965 
-0.9071 
-0.9110 
-0.9150 


0.0376 
0.0188 
0.01252 
0.00627 


90.0±0.3 
90.2±0.6 
90.6±0.8 
90.6±1.6 


6.0 


0.093 


20 a x 48 


600 


-0.8110 
-0.8170 
-0.8195 
-0.8210 


0.030 
0.015 
0.010 
0.005 


90.4±0.4 
91.0±0.7 
92.5±1.0 
95.5±2.1 










-0.7835 


0.0 


(Wilson) 


6.2 


0.068 


28 a x 56 


200 


-0.7337 
-0.7367 
-0.7378 
-0.7389 


0.021649 
0.010825 
0.007216 
0.003608 


89.1±0.8 
87.3±1.8 
86.3±2.8 
86.4±4.5 



VI. RESULTS 

A. Observations on curve fitting 

To assess the ability of the fitting algorithm to obtain 
accurately good fits to our data, we present in Figure [4] a 
histogram of the X 2 / n dof for all 320 twisted mass fits at 
each lattice spacing. 9 While the plot clearly shows that 
all the fits are reasonable, falling largely in the range 
0.9 — 1.4, the fits for the finest lattice spacing, (3 — 6.2, 
while having the same shape of distribution, have their 
mean shifted upward from that of the coarser spacings 
by about 0.3. At the finer lattice spacing more states are 
resolved so one possibility would be that the fitting algo- 
rithm has not had sufficient time to find the true minima. 
However, in comparing these fits that were allowed to go 
up to 1200 generations with those required to stop at 600 
generations, no marked improvement in the histogram is 
seen. A good fit for all our data was achieved at 600 
generations so this is not the source of the discrepancy. 

The actual source of this discrepancy can be traced 
to lower statistics, as our (3 = 6.2 data involved only 
200 configurations compared to 600 at the coarser lattice 
spacings. In Ref. 0] we fit subsets of the configurations 
of our [3 = 6.0 data to see what the effect of poor statis- 
tics would be on the result. The subfit on only 200 con- 
figurations shown in that paper clearly demonstrates an 



method 23, 24] with the fit functional form and the initial param- 
eters taken from the best fit found by the evolutionary algorithm 
on the actual data. 

These fits were for each of the 20 K PC channels at four degen- 
erate quark masses with correlators both neutral and charged, 
smeared and unsmeared. Later we also did fits to the combined E 
and T2 channels as well as for the A2 channel combining smeared 
and unsmeared correlators. There, fit quality statistics were in 
line with these, with /3 = 6.2 once again systematically higher. 
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FIG. 4: Shown are normalized histograms of the x l n dof 
(thick lines) of all the channels fit for each of our lattice spac- 
ings. To the left of the diagram one also has histograms of 
the standard deviation of these values, °' x 2/„ dof (thin lines). 
Smeared and unsmeared data at the same lattice spacing had 
essentially the same histograms so the plot does not distin- 
guish them. 



increase in the X 2 / n dof comparable to what we find here 
with our (3 = 6.2 data. Since our poorer (3 = 6.2 fits can 
be explained by fewer statistics, we conclude that where 
good fits to the data exist they are found by the fitting 
algorithm. It is notable that the shape of this histogram 
can be used to assess the minimum number of configura- 
tions required for a proper simulation. Presumably, once 
a sufficient number of configurations have been used that 
states are reliably identified, one will get a good distri- 
bution of x 2 / n dof centered on one. 10 A greater number 



Obviously the number of correlators fit in each channel varied 
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of configurations would then serve to lower the error on 
the states found and narrow the overall distribution. 11 

Since the errors of the parameters were obtained by 
bootstrapping the functional form found by the evolu- 
tionary algorithm, it was also possible to calculate the 
variation in the x 2 / n dof for each fit. Figure [4] includes 
histograms of the standard deviation a x i /„ do/ of our fits. 
For the coarser lattice spacings with 600 configurations 
one finds that there is a relatively narrow distribution of 
the x 2 l n dof across bootstrap configurations. This gives 
some confidence that the fitting algorithm has found a 
stable functional form for the fit. For (3 = 6.2, how- 
ever, the variation is found to be wider. That this is 
due to the fewer configurations available at this spacing 
can also be confirmed by consulting Ref. Q which shows 
that our subfit to 200 configurations displayed a simi- 
lar trend. Similar findings are obtained by plotting the 
standard quality of fit Q, as discussed in Ref. (3l| . 

Because our evolutionary algorithm fitting function 
was designed to fit correlators that are sums of decay- 
ing exponentials only, this causes a systematic error in 
the Af + channel where, due to a quenching artifact, this 
is not the functional form of the correlator. In quenched 
lattice QCD one has a ghost contribution in the scalar 
correlator due to the rj'-n intermediate state being light 
and of negative norm in this approximation [32| . In our 
data this effect is most pronounced at the lightest quark 
mass in the (3 = 6.2 charged channel. See Figure [5] where 
our five correlators in this channel are plotted. Because 
we only have good C and product PG in our neutral and 
charged channels respectively, this artifact will appear in 
other A% channels [25|]. Because the coefficients in our fit 
functions are constrained to be positive, the effect of neg- 
ative ghost contributions is primarily to produce a poor 
fit. Several of the higher outliers of x 2 / n dof in Figured] 
are readily traced to A\ channels containing this ghost 
contribution. 

In order to evaluate the consistency of the fitting algo- 
rithm we did at least two runs in each channel and found 
no significant variation in the results. Where difference 
occurred it would tend to be a single state with large er- 
ror found in a single correlator. Given a finite number of 
generations in a run this is to be expected. 

The fitting algorithm appears to have resolved states 
predictably. Unsmeared channels found greater numbers 
of states than did smeared channels. In some channels 
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greatly so the width of the distribution is partly related to the 
statistics induced by this variation. 
11 One may wonder why our smeared and unsmeared correlators 
had comparable histograms if the smeared data had to resolve a 
fewer number of states. Why this should be the case is likely due 
to the asymmetrical smearing that is done to our source and sink 
operators. While it is true that we are removing excited states by 
smearing, in our case we are producing noisier correlators with 
the result that our statistics remain wanting for our smeared 
correlators at our finest lattice spacing. Our baryon study with 
the same smearing showed a similar result [2{| . 




FIG. 5: Shown are the five unsmeared diagonal correlators 
for the charged Af + channel at f3 = 6.2 and lightest quark 
mass. Solid points are positive values and open are negative. 
In parentheses after each number the octahedral irrep of the 
local quark structure and the gauge field used to generate the 
operator are given. The first operator has no extended gauge 
field structure. The last two operators use different T" local 
quark representations in their construction, operator #4 using 
quark structure T-j 2 , with the superscript 2 indicating the 
multiplicity a from Table IVl 



no state was found, which affirms that the method does 
not claim signals in everything. Occasionally a clearly 
spurious state with large error would appear in a sin- 
gle correlator which is also expected from the statistical 
variation of the data itself. Finer lattice spacings in the 
unsmeared case found more states since a greater num- 
ber of excitations could be resolved. Occasionally the 
[3 = 6.2 lattice spacing appeared to resolve an interme- 
diate state not seen in the other spacings. While the 
bifurcation of states could certainly occur at finer spac- 
ing, the fewer configurations at j3 = 6.2 do not allow us 
to make a strong conclusion from this observation. 

Occasionally states of low energy would appear in some 
channels with large error below the obvious ground state. 
Such errant signals are a consequence of fitting at Eu- 
clidean times far from the source, where correlators are 
noisy. They may be an artifact of time step correlation 
not being entirely removed through the use of the corre- 
lated x 2 ■ We regard the presence of these noisy spuri- 
ous states as a reminder that we are using a true black 
box method: the X 2 / n dof is slightly reduced by their 
presence, and we have not prevented their appearance 
through human intervention such as choosing a fitting 
window in Euclidean time. 

As found in Ref. [6], the computer time required for 
an evolutionary fit scales with the number of parameters 
required. This depends not only on the number of states 
in the channel but includes the coefficients required on 
each correlator, and hence depends on the number of 
correlators fit in the channel. By smearing we reduce the 
number of states and therefore the number of coefficients 
and thus smearing allowed us to fit many correlators in 
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the same channel quickly. 

As part of the analysis we tried to fit all smeared corre- 
lators corresponding to the J = 4 channel which includes 
all four octahedral irreps {At, E, T%, and T2) to which its 
subduced representation reduces for a given PC. Such 
fits involved fitting on the order of forty correlators simul- 
taneously and, although taking proportionately longer, 
were, in the end, tractable by the fitting algorithm. 

B. Observations on twisted mass 

Figures [5] and [JJ show fits with both charged and neu- 
tral operators in the A% and T\ channels, which couple to 
spin zero and one respectively in the continuum. Results 
at four different degenerate quark masses are shown. The 
sizes of the points are scaled to give a qualitative impres- 
sion of the significance of the state in the datasets via 
the factor: 




Here N is the number of correlators involved in the fit, 
is the coefficient on the i th correlator for the given 
energy state (taken to be zero if there was no coefficient 
found on that dataset) and e'*' is the bootstrap error of 
the given coefficient. The arctan function is used to en- 
sure the scale factor ranges from to 1. The utility of 
the scaling is that it facilitates the identification of the 
same state across lattice spacings and quark masses for 
extrapolation purposes. The scale factor is not numeri- 
cally involved in any extrapolation however. 

In Figures [6] and [7J one is able to see easily the effect 
of the twisted mass channel contamination discussed in 
Section IIIII In Figure [6] the ground state pseudoscalar 
("pion") is clearly visible across quark masses where it 
should be in the + channel for both the charged and 
neutral correlators. For the charged operators, channels 
are expected to contaminate others with the same PG 
product. The charged pion shows up clearly as con- 
tamination in the A^ channel as predicted. That it 
is contamination is clear not only because there is no 
corresponding signal in the neutral channel but also the 
larger error and smaller symbols are both indicative of a 
weaker signal than the authentic one in the A\ + chan- 
nel. Comparison of the A^ + channel with the A\ chan- 
nel similarly shows contamination of the authentic scalar 
A^ + ground state signal in its twinned channel. 

Turning to the neutral case, charge conjugation is re- 
spected and one sees the pion now contaminating the 
A~l + channel. States in the latter channel in turn are 
contaminating the A^ + channel where one clearly sees 
weak contamination states interspersed with authentic 
states. It is a testament to the fitting algorithm that it 
is able to distinguish these errant states. 

The strength of the ground state vector ("p meson") 
allows us to confirm similarly the contamination relations 



for those channels not evident in the scalar case. The p is 
clearly identified in both charged and neutral T x chan- 
nels in Figure [JJ Consideration of the contamination in 
the charged channel suggests the particle should appear 
in the channel, which it clearly does. In the neutral 
case the p contamination is appearing in the channel 
as expected. It is worth noting that contamination states 
are not appearing in unexpected channels. The strength 
of the pion and the p as signals, even in contamination, 
would make them easily discernible. 

When contamination is weak or when the actual states 
in the channel being contaminated are themselves weak 
one may see distortions in the actual state if the two 
states are nearby or cannot be distinguished. Higher 
statistics are required to disentangle correctly the greater 
number of states that will appear in a given channel. The 
value of being able to compare neutral and charged chan- 
nels which get contaminated differently is readily appar- 
ent. 

Overall, contamination appears to decrease with finer 
lattice spacing. At lighter quark mass the contamination 
also lessens but so does the actual signal in the original 
channel. In our simulation, smearing reduced the con- 
tamination strength and in our smeared fits it was only 
possible to identify reliably ground state p and pion con- 
tamination. The decrease in contamination with lattice 
spacing and the effect of smearing can be seen in the fit 
plots in Section IVI CI The fits in Figures [5] and [7J be- 
ing to data both unsmeared and at our coarsest lattice 
spacing, are for illustration of the contamination effect. 
The smeared data we actually fit for our results appear 
to have limited contamination identified by the fitting 
algorithm. 



C. Meson spectroscopy 

The data from the smeared T" and T^~ + channels 
are provided in Figures [8] and [9] respectively. The sizes of 
all data points are scaled by Eq. ([55)1 so tiny points are 
often devoid of any physics content. Similar plots for all 
other A PC appear in Ref. (3l| . 

In Figure the extrapolated ground state vector me- 
son (the p) is in agreement with its physical value, and a 
single excited state is resolved. The neutral (3 = 5.85 
excited state, while present at every quark mass, ap- 
pears to be interfering with the ground state. The re- 
sulting poor chiral extrapolations of the neutral (3 — 5.85 
data will be excluded from the continuum extrapolation 
of both the ground and excited states. At the heavi- 
est quark mass the ground state for the neutral (3 = 6.0 
fit is bifurcated (the points are indistinguishable on this 
plot) presumably because statistical variation made this 
a slightly more probable fit. Due to limited statistics in 
the /? = 6.2 channel (recall Section fVI A|) . there may be 
increased systematic errors that arise during the fitting 
at (3 = 6.2. One might therefore prefer to omit (3 = 6.2 
data from continuum extrapolations. 
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FIG. 6: The results of the fits to the four A\ channels at j3 = 5.85 for the unsmeared correlators. For each of the charged (C) 
and neutral (N) channels, there are four fits for the quark masses going from the heaviest (mass #1) to the lightest (mass #4). 
Points are scaled to the p mass found in the smeared charged mass #2 channel. Point sizes are scaled to reflect their significance 
as given by Equation (|35[) . 
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FIG. 7: The results of the fits to the four Ti channels at /3 = 5.85 for the unsmeared correlators. Notation is identical with 
Figure \§\ 
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FIG. 8: Fits to smeared T 1 correlators are shown. Experimental measurements [331 ] are plotted for comparison. 
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FIG. 9: Fits to smeared T' 1 _+ correlators are shown. A clear signal appears between 2 and 3 GeV which has been extrapolated. 
The lower points seem to suggest that the fitting algorithm is detecting some other state, likely channel contamination due to 
tmLQCD. This is an exotic channel. Experimental claims [33J are plotted for comparison. 
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FIG. 10: Experimental values |33( | of light unflavored mesons 
(squares) in the isovector channel for each J PC are com- 
pared to the results of this simulation (diamonds) that are the 
weighted average over lattice spacings /3 = 5.85 and f3 = 6.0 
extrapolated to zero quark mass. 



In Figure [SI results in the exotic Ty channel are plot- 
ted. A clear signal is evident between 2 and 3 GeV, while 
tiny points also appear in the small mass range. These 
tiny points may be attributed in part to contamination 
with the other T\ channels due to tmLQCD's symmetry 
structure, and in part to statistical variation within the 
fitting procedure. At our largest quark mass, which is 
approximately the strange quark mass, our results are 
consistent with the lattice study of Ref. @ which used 
displaced quarks to generate angular momenta in the op- 
erators for quenched simulations at (3 = 6.0. 

The results of all our simulations, after linear chiral 
extrapolations, are given in Table IXIVI along with con- 
tinuum extrapolations. Extrapolations that omit the 
(3 = 6.2 data, i.e. simple averages of data at the other 
two f3 values, are provided in Table IXVI for all A PC . 
The predictions for charged meson masses are also shown 
graphically in Figure 1101 together with the experimental 
mass spectrum. 

Note that our results from the heaviest quark mass at 
(3 = 6.0 (provided explicitly in Ref. [3lj |) can be com- 
pared directly to Ref. 0], and good agreement is found 
for the available channels: ao, 7T, p, &i, a± and a,2- The 
exotic states which we identify in the H , 1 ^ and 2" 1 
channels, for the strange quark mass at (3 — 6.0, are also 
found to be comparable to values which they ascribe to 
those same J PC . 

Our result for the 1 h from Table IXVI agrees well with 
the general consensus of quenched [9l.ll6l [34 . HH, [36| and 
dynamical [33, [38| lattice studies that the lightest exotic 
meson is 1 h with a mass near 1.9 GeV. We also note 
that the flux-tube model predicts exotic hybrids H , 
1 h , and 2" 1 near 1.9 GeV in the light isovector chan- 
nel [39j |. 

Although results of our simulations with Wilson 
fermions have not been reported here, they did provide 



independent confirmation of which signals were actually 
contamination that arose due to the twisted mass term. 

In all but the A2 channels, the smeared operators pro- 
duced better fits than unsmeared operators. For A PC , it 
was found that just fitting the smeared channel tended 
to produce failures in the bootstrap fits which are used 
to generate the errors of the parameters. The reason 
for this failure is poor statistics due to a combination of 
the relatively few correlators in each A2 channel (no row 
averaging occurs in a one-dimensional channel) and the 
weakness of the expected states in the channel since A2 
couples only to continuum J > 3 as shown in Table [II] 
As our motivation for fitting smeared states is to reduce 
the influence of high energy states and since the A2 chan- 
nels have very few low lying states even in the unsmeared 
case, we refit all A2 correlators, both smeared and not, 
simultaneously. This gives sufficient statistics to make 
our bootstrap fits converge properly. 

Visual inspection of the fits to the E and T2 channels 
independently, in each PC combination, revealed that re- 
sults were remarkably similar for any state found. This 
is to be expected since, by Table [Til the lowest angular 
momentum to which they both couple is J — 2. With 
the assumption that the few states found across these 
channels arc of this angular momentum, we do combined 
fits of all correlators of E and T2 together for each PC. 
In these combined fits the algorithm finds the states to 
be statistically the same (i.e. it does not bifurcate any 
of the states), thus allowing us to conclude that our re- 
sults are commensurate with ascribing a value J = 2 to 
the common state. 12 Furthermore it allows us to use all 
the available data to extract its value. We thus tabulate 
results only for the combined (E, T%) fits, and not them 
separately. Our lattice spacing is fine enough that given 
our statistics it is unlikely that any splitting between E 
and T2 states for a common J = 2 state would be de- 
tectable, and it is not. 13 The remaining states in our 
fits are assigned their most probable continuum J values 
based on Table HTl namely J = 0, 1, and 3 for irreps A\, 
Ti, and A 2 respectively. 14 



D. Observations on operators 

To give insight into the nature of the operators which 
contribute to a particular state we have produced Ta- 
ble IXVII in Appendix [A] There one finds the operators 



12 Technically J = 4 is also a logical possibility since its subduccd 
representation reduces to both irreps E and T2 as well. 

13 Such splitting can be observed across continuum subduction 
channels on coarser lattices. See Ref. |4dl |. 

14 While J = 3 could also appear in T\ and T2 channels as well, the 
likely weakness of such a signal made a combined fitting of these 
channels inappropriate in comparison to the J = 2 case, given 
the presence of lower J in Ti and T2 channels and the limitations 
of our statistics. 
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TABLE XIV: Extrapolations of each resolved channel labeled by its irrep (A ) and its type, charged (C) or neutral (N), are 
given. States designated by an asterisk indicate a higher energy state in the same channel. The inferred continuum quantum 
number J of the state is also given. Chiral extrapolations (linear) for each /3 are shown along with their x 2 / n dof in subsequent 
columns. The continuum extrapolation of these results is shown in the last column. We have omitted the pion (A± + ) from the 
table. All energies are in GeV. 



Channel 
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5.85 


P = 


6.0 
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6.2 


Final 


Irrep (s) 
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C/N 


Energy 


X 2 /v 


Energy 


X 2 /v 


Energy 




Energy 


X 2 /v 







C 


1.35(3) 


1.579 


1.36(3) 


0.131 


1.40(6) 


0.015 


1.40(6) 


0.231 






N 


1.01(4) 


2.566 


1.21(4) 


0.198 


1.20(8) 


0.012 


1.39(8) 


1.683 







C 


3.3(3) 


0.015 


4.2(4) 


1.576 


5.1(10) 


0.023 


5.6(8) 


0.053 






N 


3.5(4) 


0.688 


2.5(3) 


0.014 


3.1(6) 


0.194 


2.1(7) 


2.209 


i 





C 


0.31(17) 


2.815 


1.96(18) 


0.048 


2.08(10) 


0.303 


2.2(3)° 








N 


1.47(16) 


0.027 


1.69(10) 


0.268 


1.96(13) 


1.548 


2.14(19) 


0.309 


T++ 


1 


C 


0.98(5) 


0.371 


1.25(5) 


0.990 


1.6(2) 


0.376 


1.65(12) 


0.461 






N 


1.33(2) 


0.101 


1.428(17) 


0.003 


1.52(10) 


0.010 


1.57(5) 


0.061 


T*++ 


1 


C 


2.4(3) 


0.843 


2.3(3) 


0.213 


3.3(6) 


0.315 


2.9(6) 


1.688 






N 


3.0(7) 


0.109 


2.4(3) 


0.298 


5.8(13) 


0.459 


2.9(10) 


7.283 


T+- 
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1.40(3) 


1.507 


1.48(3) 


1.018 


1.74(8) 


0.961 


1.69(7) 


4.136 
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1.14(6) 


0.715 


1.62(6) 


0.018 


1.57(7) 


0.199 


1.90(10) 


6.728 


T *+- 
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2.6(4) 


1.095 


2.7(3) 


0.535 


2.8(6) 


0.011 


2 9(7) 


0.007 






N 


3.1(4) 


0.283 


2.8(4) 


0.021 


3.2(4) 


1.833 


3 1(6) 


0.591 


± i 
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1.80(18) 


0.774 


2.15(18) 


0.008 


2.4(4) 


0.258 




0.022 






N 


1.7(2) 


0.566 


2.3(2) 


0.027 


2.9(5) 


0.000 


3.3(5) 


0.157 


T — 


1 


C 


0.767(19) 


0.141 


0.764(16) 


0.036 


0.79(3) 


0.066 


0.78(3) 


0.573 
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0.797(6) 


0.508 
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0.290 


1.0(4) 


0.244 


1.3(4) 


0.003 






N 


0.6(4) 


0.049 


0.52(13) 


0.179 


0.5(3) 


0.051 


0.4(4) 


0.002 


(E,T 2 )*++ 


2 


c ' 

L> 


3.7(7) 


ft i of; 


2.7(4) 


U.Uoo 


2.y(y J 


U.lz ( 


1.9(10) 


0.468 






N 


2.4(3) 


2.235 


3.0(4) 


0.103 


3.2(8) 


0.037 


3.7(8) 


0.023 


(E,T 2 )+- 


2 


C 


2.8(6) 


0.057 


2.7(4) 


0.065 


2.9(8) 


0.018 


2.8(9) 


0.066 






N 


3.2(5) 


0.005 


2.9(4) 


0.031 


3.1(6) 


0.026 


2.9(8) 


0.172 


(E,T 2 )~+ 


2 


C 


2.3(2) 


0.334 


2.9(4) 


0.025 


3.3(9) 


0.121 


3.7(7) 


0.001 






N 


2.4(2) 


4.133 


2.8(4) 


0.001 


3.0(7) 


0.344 


3.2(7) 


0.001 


Ai + 


3 


C 


0.5(7) 




1.8(5) 


0.002 


3.5(8) 


0.052 


4.3(10) 


0.573 






N 


0.7(4) 




10.(8) 


0.005 


3.2(8) 


0.017 


4.4(11) 


0.912 




3 


C 


1.3(4) 


0.230 


1.4(7) 


0.672 


1.1(6) 


0.072 


1.1(8) 


0.113 






N 


1.5(5) 


0.376 


1.2(4) 


0.088 


1.2(4) 


0.664 


1.0(6) 


0.009 



"Continuum extrapolation excluded (3 = 5.85 data point. 



which are found to contribute to each extrapolated state, 
written in terms of their local quark and extended gauge 
field content from which they are constructed. The na- 
ture of our operators is such that only one irrep is found 
in any such product as shown in Table |X] so this labeling 
uniquely identifies the operator. 

To quantify the significance of an operator's diagonal 
correlator to a given channel two additional pieces of in- 
formation are given in Table IXVll For each such oper- 
ator the number of times it was found contributing to 
the state out of the twelve possible lattice spacing and 
quark mass combinations that were fit is given. In some 
cases a small contribution may not be statistically sig- 
nificant and as such the fitting algorithm may opt not 



to give a coefficient for it. Secondly, the table provides 
the largest magnitude of the coefficient Z found for that 
state, among those fits that had a coefficient. These lat- 
ter values are scaled by their error so as to give a measure 
of their significance. The table is organized to give the 
most significant operators contributing to a channel first 
in the list. 

Several observations may be made from Table IXVII 
For one, operators with a non-scalar gauge structure 
still contribute to ground states. For the p meson (Tj" 
ground state) the greatest contributions come from those 
operators which have a quark structure of , either 
a local operator or with an extended gauge field. 

However, as can be seen in the table, one also has a signif- 
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TABLE XV: Weighted average of the chirally extrapolated 
[3 = 5.85 and j3 — 6.0 lattice spacing results of each resolved 
channel labeled by its irrep (A PC ) and its type (charged or 
neutral) are given. States designated by an asterisk indicate a 
higher energy state in the same channel. The inferred contin- 
uum quantum number J of the state is also given. We have 
omitted the pion (A± + ) from the table. All energies are in 
GeV. 



Channel Energy 



IlTGp(s) 


j 


Charged. 


Neutral 


4++ 


n 

u 




i 

1 . IZ^O j 


A+- 


n 






A*,~ + 
"■i 


o 


l q6d8i a 


1.63(8) 


T++ 


1 


1.12(3) 


1.394(14) 


T *++ 


1 


2.4(2) 


2.5(2) 


T+- 


1 


1.44(2) 


1.37(4) 


T*+- 


1 


2.7(2) 


2.9(3) 


r f + 


1 


1.98(13) 


1.98(16) 




1 


0.765(12) 


0.812(12) 


r*- 


1 


1.42(6) 


1.63(10) a 


(E,T 2 )++ 


2 


0.61(15) 


0.53(13) 


(E,T 2 )*++ 


2 


2.9(3) 


2.6(3) 


(E,T 2 )+- 


2 


2.8(3) 


3.0(3) 


(E,T 2 )-+ 


2 


2.46(19) 


2.52(20) 


Ai + 


3 


1.3(4) 


0.7(4) 


A 2 + 


3 


1.3(4) 


1.3(3) 



"Excluded = 5.85 data point. 



icant overlap with states of the form A± ®T± in which 
the gauge field is actually providing the vector nature of 
the state. 15 The excited state T* shows, on the other 
hand, a much greater relative contribution from the vec- 
tor gauge field operators, but still dominant contributions 
from the vector quark operators. For the ground states 
T 1 ++ and T± the operators with quark vector struc- 
ture dominate, however for their excited T 1 * ++ and 
states it is seen that a greater contribution comes from 
the operators with vector gauge field structure. This ob- 
servation is an example of how lattice studies of the spec- 
trum can inform discussion of the distinction between hy- 
brid vector mesons and conventional vector mesons which 
is made in models like the flux-tube model [4l| . 

Table IXVII also grants further insight into the geome- 
try of related states. In the constituent quark model, for 
instance, the an, ai, and b\ states belong to a common P- 
wave and hence are considered to share common angular 
momentum [9]. If one looks at the main contribution to 
these states in their corresponding octahedral channels, 
A^ + , T^~ + , and T 1 _+ respectively, it is observed that af- 
ter the local operators the main contributor in each of 



these channels is of the form X PC ® T± , where X PC 
is whatever local spin structure is required to produce 
the correct channel. This suggests that just as common 
spatial angular momenta (L = 1) are traditionally con- 
sidered to unite these channels, so too does similar gauge 
field structure. 

For those states we identified as J = 2 with our com- 
mon fits to E and T 2 channels, Table IXVTl shows further 
affirmation of this identification. In the continuum, five 
J = 2 operators can be formed with the coupling of a 
,7=1 quark and J = 1 gauge field structure. On the lat- 
tice these five states still appear in a coupling of T± quark 
and Ti gauge field, but now in the subduction of J = 2 
as a two-dimensional E irrep and a three-dimensional T 2 
irrep. As such one expects these E and T 2 operators 
produced via T\ ® Ti corresponding to an actual J = 2 
state to have similar properties. 16 A consideration of our 
eight (E, T 2 ) common fits shows that the contributions 
of such E operators (denoted by T\ ® Ti) have compara- 
ble counts and significance to their related T 2 operators 
(denoted by T\®Ti) in all cases. Only in the compa- 
rably weak 2 ++ ground state fit, and there only in the 
neutral channel, is the pattern less obvious. For all these 
states the importance of the gluonic degrees of freedom is 
readily apparent, in contradistinction to the limitations 
imposed by a simple quark model for two local quarks. 



VII. CONCLUSION 

Simulations with twisted mass lattice QCD (tmLQCD) 
have been used to explore the spectrum of mesons hav- 
ing all possible values of angular momentum, spin, and 
parity. Numerical results are consistent with those ob- 
tained by other authors using other lattice actions. As 
seen in Figure 1101 we are still far from seeing the number 
of mesons claimed experimentally, and uncertainties are 
still large, but a viable methodology has been put into 
place. In particular, questions related to the appropriate- 
ness of tmLQCD for this physics can now be answered, 
the usefulness of operators containing local quark and 
anti-quark fields can be evaluated, and the valuable qual- 
ities of an evolutionary fitting method can be confirmed. 

Because tmLQCD breaks the parity symmetry of the 
strong interactions, it is not clear a priori whether it can 
be used to determine the spectrum of hadrons having def- 
inite quantum numbers J PC , which appear as A PC on a 
discrete lattice. In the case of both isovector and isoscalar 
mesons, we have identified orthogonality relations among 
certain operators which are valid in tmLQCD for any 
chosen twist angle. We have also emphasized that some 
states cannot be separated, such as A ++ ,A in the 



5 For the ground state pion not shown in Table IXVTl 

it is 

to be remarked that there is similarly a clear ground state signal 



from the T, 



operators. 



16 However, since one is row averaging over three operators in the 
T2 case compared to two in the E case the correspondence is not 
expected to be identical due to differing statistics. 
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isovector case and A ++ , A h in the isoscalar case. Nev- 
ertheless, simultaneous fits to multiple correlators allow 
the complete spectrum to be obtained (as proposed in 
Ref. j4j), and comparison of the isovector spectrum with 
its isoscalar counterpart serves to clarify the PC quan- 
tum numbers of the states. During this work, simulations 
with the Wilson action were performed also in order to 
confirm these claims. 

Meson operators that require only one quark prop- 
agator, i.e. operators with a common local source for 
mass-degenerate quark and anti-quark, minimize compu- 
tational expense. Connecting the quark and anti-quark 
to spatially extended gauge fields permitted all possible 
meson quantum numbers, but how well do these opera- 
tors overlap with the physical states of QCD? Our numer- 
ical explorations produced signals for both conventional 
and exotic mesons. The strength of each operator's over- 
lap was tabulated, so that future studies can make in- 
formed choices of operators for the various meson chan- 
nels. Gauge field smearing was included at source and 
sink, while quark field smearing was performed only at 
the sink. 

The use of an evolutionary fitting algorithm avoided 
the danger of human bias during data analysis. The algo- 
rithm was able to fit multiple correlators simultaneously, 
with some states shared (or not) across data sets. The al- 
gorithm itself determined how many states were present 
based on statistical signficance, and all time steps beyond 
the source were included in all fits. Survival of the fittest 
was defined to mean minimizing X 2 / n dof- The algorithm 
even served to alert us to a case of insufficient statistics 
in one ensemble of lattices; this information appeared as 
larger x 2 / n dof values output from the fit, which could 
only be reduced toward unity by increased statistics. 

Based on this work, we conclude that a detailed study 



of the full meson spectrum using tmLQCD is feasible. 
However, to obtain precise results, substantially more 
configurations would be required than we have used, 
and additional classes of operators should be considered. 
Quark smearing at the source, to match the sink, should 
improve the signal-to-noise ratio. With more precise me- 
son mass determinations at a few f3 values, one could 
see directly whether the continuum extrapolation of tm- 
LQCD offers a significant advantage over other lattice 
actions, commensurate with the drawback of tmLQCD's 
lack of parity conservation. It would also be interesting 
to consider the viability of using tmLQCD for a thorough 
study of the baryon spectrum. 



APPENDIX A: OPERATOR CONTRIBUTION 

The operators contributing to each extrapolated state 
are listed in Table IXVIl Operators are identified 
by the direct product of quark local irrep and ex- 
tended gauge field irrep from which they are projected, 
pat f i ^ \P u c u _ Operators without an extended gauge 
structure given are local quark operators. Quark irreps 
from the second multiplicity have a superscript 2. In the 
case of the (E,Ti>) fits we designate T2 operators with 
an <g> and E operators with (g). In the A 2 fits where 
both smeared and unsmeared correlators were fit the un- 
smeared correlators are distinguished similarly with ®. 
After each operator the number of lattice spacing-quark 
mass channels in which a coefficient was identified for the 
state is given (maximum 12) followed by the maximum 
value of the coefficient Z divided by its error over those 
channels. 



TABLE XVI: Table of operators contributing to each extrapolated state. 



Channel 




Irrep(s) C/N 

A++ n A++tm rvc on A++ < 


Contributing Operators 

> A++tin rM a\ *>.rn+-m o m'2 ^ rr, /oo^oN 



N A ++ ® A++(12, 48.4), A++(12,48.1), T+" ® T+-(9,4.34), Tf— ® Tf- (3, 3.92) 
A+~ C A+"(12,9.25), A+~ ® A++(12, 9.25), Tf+ <g> T+ (10, 5.74), Tf~~ ® Tf +(7,4.60), 

T~ ®Tf+ (4, 4.41) 

N A+-(12,5.26), A+~ ® A++(12, 5.20), Tf+ ® Tf" (11, 4.19), Tf ~ ® Tf+ (5, 3.11), 
Tf~~ ®Tf+(2,1.30) 

A{-+ C Af+(8,29.5), A-+ ® A++(8,28.9), Tf" ® Tf"(8, 11.3), Tf— ® T+~(8, 3.57), 

T+~ ®T~ (4,9.21), Af-+(4,3.11), T++ ® Tf+(4, 0.92), A\~+ ® A++(3, 3.46) 
N A± + ® A++(12,26.9), A^+(12, 26.9), A{~ + ® A ++ (10, 5.53), Tf ® T +_ (8, 3.26), 
A?-+(7,5.58), Tf" ® Tf" (6,4.38), Tf" ® Tf"(5, 5.08), Tf+ <g> Tf +(4,4.92) 
T++ C T++® A++(ll,36.8), Tf +(11, 36.7), Tf" ® T+"(7, 2.47), T+~ ® E+'(7, 2.00), 

Tf" ®T 2 — (3,2.14), Tf~~ ® Tf" (3,2.12), A+~ ® Tf" (3, 1.32), Tf~~ ®T x (2,1.45), 
Tf" <g> Tf" (2, 1.39), T++ ®T++(2,0.60), T++ ® E ++ {1, 1.42) 
N T++ ® A++(12,61.8), T++(12, 50.0), T+" ® T+-(8, 5.43), Tf" ® E+~ (7, 3.54), 

If— ®Tf- (7, 1.69), A+- ®T+-(6, 1.97), Tf ® T 2 (5, 4.16), T{~ ® Tf" (5, 2.98), 
Tf" ®Tf- (5,0.80), T++ ® £++(2,0.69), T++ <g> T++(l, 0.77), A~+ ® Tf+(1, 0.03) 
T*++ C A+" ®Tf _ (12, 6.56), Tf" ®T+-(11, 6.70), T++ <g> A++(8, 9.12), T++(7, 9.46), 

Tf- ®Tf- (6,5.24), Af+ ® Tf +(6, 2.76), Tf— ® Tf- (6,2.55), T~ ® Tf" (6, 1.96), 

Af- + ®rf + (5,4.11), T++ ®T++(5,2.46), T++ g £++(5, 1.68), Tf~~ g Tf" (3, 4.45), 

Continued on next page 
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Channel 

Irrep(s) C/N Contributing Operators 

,3.53), 

>T 2 — (7,1.31), 
Tf" (4,2.47), 

69 £/ (4, 1.05) 

" <8> A++(12,38.6), T+-(12,35.3), T++ <8> T+"(8, 6.19), A++ ® T+"(7, 2.78), 
A7/+ ®Tf "(5, 2.56), T++ <8> £+"(5, 2.41), T x ® r 2 ~ + (4, 3.83), Tf— ® Tf+(4, 1.59), 
^-+(g)T 1 ~(4,1.41), Tf"" ®T 2 -+(2,1.64), T+" ® T++(2, 0.98) 
N Tf _ <8> ^++(12,29.9), Tf-(12,29.6), Tf+ <8> Tf _ (8, 4.67), Af+ ® Tf _ (8, 3.92), 

Tf" <8>Tf+ (3, 3.44), ®Tf- (3, 2.34), T++ <8> £+"(3, 2.30), Tf" <g> Tf+(3, 1.16), 

Tf _ <8> T 2 ++ (2, 1.16), (8) Tf _ (2,0.97), Tf" ® T 2 "+(l, 0.96), Tf ® Tf+(1, 0.67) 
Tf+ _ C T++ ®Tf-(10,7.08), A++ <g>Tf-(9,5.89), Tf - ® A++(9, 4.79), Tf _ (9, 3.68), 

Tf" ® Tf +(7, 2.83), <8> Tf" (6,4.69), A"+ <g> Tf "(5,4.93), T++ <8> £+"(4, 2.22), 

Tf— <8>Tf+(3,3.24), Tj ® Tf+(3, 2.48), <8> Tf +(3, 2.17), Tf" ® £++(2, 2.88), 

T+- <8> T++(2,0.88) 

N T+ + (g>T 1 + -(10,7.28), A++ (8)T+-(9,7.23), T+"(6, 9.53), T+" ® A+ + (6, 9.52), 

<8> Tf" (5, 5.36), Tf" <8>Tf+(5, 5.30), A~+ ® Tf" (4,4.53), Tf -- ® Tf+(4, 4.02), 
Tf" ® T 2 -+(4,3.27), Tf— ® Tf + (4, 2.82), T++ <8> £+"(4, 2.59), T+" <g> £++(3, 3.36), 
Tf" <8>T 2 ++(3, 1.73) 

Tf+ C Tf"" (8)T+-(12,7.54), Tf" ® Tf" (12, 7.31), T+" ® T 2 ~(4, 3.99), Tf~ ® E+~(4, 2.73), 

Tf" <8>£+ _ (4, 2.56), A+- (g>Tf-(4,1.66), T++ <8> Tf+ (3, 3.67), A++ ® Tf +(3, 1.81), 
Tf " <g> Tf "(3, 1.36), Tf + ® Tf +(2, 2.70) 
N Tf" ® Tf" (11, 7.29), Tf -- <g> T+-(ll, 7.09), T+" ® Tf "(9, 3.45), T++ ® Tf +(7,3.03), 
A++ ® Tf +(6,2.78), A+- ® Tf" (4, 4.50), T++ ® Tf+(4, 3.87), Tf ~ <8> £+" (4, 2.34), 
Tf - ®£+- (4, 1.97), T+- ®Tf- (3, 3.02) 
Tf- C Tf-- <g) A++(12,64.1), Tf— (12,64.1), Tf" (12, 43.7), Tf" <g> A++(12, 43.7), 

® T+~ (7, 7.34), A7/+ ® Tf "(6, 4.84), Tf ~ ® E++(4, 2.43), Tf" ® T++ (4, 0.95), 
T++ ® Tf- (3, 1.94), T++ ®T 2 ~ (3, 1.85), Tf~ ®T++(3, 1.50), A++ ® Tf" (3, 1.29), 



N 



Tf - ® £++(3,0.84), Tf - ®Tf+(2, 2.74), A+~ ® Tf +(1, 1.32) 
Tf- (12, 107.), Tf - ®A++(12, 107.), Tf-- (10, 25.7), Tf~ <g) ^(10, 25.1), 
A7/+ ® Tf "(9, 8.06), A\~+ ®T+~(8,7.07), T++ (g) Tf" (8, 1.35), A++ (g> Tf" (6, 3.08), 
Tf - ® Tf +(5, 2.42), T++ <g> Tf" (5, 2.41), Tf " ® £++(5, 1.22), Tf ~ <8> £++(3, 1.87), 
A+- ® Tf+(3, 0.41), Tf - ® T++(2, 1.86), T+" ® Tf+(2, 1.83), Tf-- <8> T++(2, 1.00) 
Tf"" C Tf- (12, 22.0), Tf - <g> A++(12, 21.9), Tf "(11, 17.4), Tf-- ® A++(ll, 17.3), 

A^+ ®T+-(11,4.39), (g>T+-(10,7.17), T++ ® Tf "(7, 2.83), A++ ® Tf "(7, 2.81), 

Tf - <g> £++(7,1.63), T++®T 2 — (5,2.03), Tf" ® T++(4, 1.05), Tf" ® Tf +(2, 1.57), 
T+- <8>T 2 -+(2, 1.43), A+- ® Tf +(2, 0.60), Tf ~ ® £++(1, 1.62), Tf — <g> T++(l, 1.01) 
N Tf-- ® ^++(12,46.9), Tf— (12,45.3), A 2 -+ ® Tf" (10, 4.82), Tf - ® A++(9, 19.3), 
Tf- (9, 18.9), A++®Tf- (7,3.16), Af+ <g> T+-(6, 4.15), T++ ® Tf" (6, 4.07), 
Tf-- ® £++(3, 1.48), T++ <g)T 2 — (3,1.30), Tf" ® Tf+(3, 1.25), A+" ® Tf +(3, 1.02), 
Tf-- <g>T 2 ++(l, 1.72) 

(£,T 2 )++ C T++®£++(7, 1.09), Tf ~ jg> T 2 ~ (6, 1.80), Tf~ <g> Tf" (6, 1.24), Tf" ®T 2 ~ (6, 1.13), 
Tf~(g)Tf-(4,1.79), T+-®T+- (4, 1.24), Tf " <g) T 2 ~ (4, 1.07), A+" ® £+" (4, 0.99), 
T+- (g)T+-(3, 1.39), A- +®Tf+(3,0.99), T+-®A+-(3, 0.88), Tf - ® Tf- (2, 0.74), 
A++®T++( 1,1.16), Tf-^Tf- (1,0.55), Tf~_®T 2 ~ (1,0.38), T+-®£+- (1,0.04) 
N T+-®£+-(6, 1.35), Tf-®Tf-(5, 1.89), T+-®T+-(4, 1.55), Tf" ®Tf (4, 1.48), 
Tf--_®T 2 ~ (4, 1.42), Tf"" ®Tf" (3, 2.20), T+" ® T+"(3, 1.28), T++®£++(3, 1.04), 
Tf" (8>^+- (3, 0.89), A+- ®£+-(2, 1.31), Tf " (8> T 2 ~ (2, 0.94), Tf ~ ® T 2 ~ (2, 0.93), 
A++®T++(2,0.69), T++ ® T++(l,0.74) 
(£,T 2 )*++ C T+-®T+- (11, 5.66), Tf" <g> T+-(ll, 4.66), A++ ® £++(7, 2.06), A?-+®T 2 -+(6, 2.27), 
Tf - ®T 2 — (6,2.06), Tf~ ® T 2 ~ (6, 2.02), Tf-®T 2 ~ (6, 1.75), A++®T++(6, 1.52), 
T++®£++(4,2.27), Tf-- ® Tf "(4, 2.23), Tf~®Tf "(4, 2.20), ^^+®T 2 -+(4, 2.20), 
Tf-®Tf- (4, 2.19), Tf" <8>Tf- (4, 2.05), Tf~®T 2 ~ (4, 1.96), Tf-®£+- (4, 1.82), 
T++®T++(4,1.71), T+-®A+-(4,1.61), T++ ® £++(4, 1.51), A+~ ® £ + "(4, 1.12) 
N T+-®T+- (12, 5.50), Tf - ® Tf" (12, 4.73), Tf "®Tf " (8, 3.08), Tf~ (g) T 2 ~ (7, 3.85), 
Tf - ® Tf- (7,3.50), Tf--®Tf- (7, 2.96), A++®T++(7, 1.92), Tf^Tf "(5, 3.85), 

A+~ (8) £+'(5, 2.91), Tf" g Tf" (4,4.34), Tf—^Tf "(4,4.23), A^+^Tf +(4, 3.69), 

Continued on next page 
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Channel 

Irrep(s) C/N Contributing Operators 

Af"+®Tf+(4,3.54), Tf-®£+~ (4, 3.03), T+~®A+~(4, 1.59), A++ ® £++(3, 3.59), 
Tf" <g> Tf" (3, 2.30), T++ <g T++(3, 1.92), r++®T++(2, 1.34), Tf+®£++(l, 2.04) 

(£,T 2 )+- C T 1 ++®T 1 + - (11, 4.69), Tf+ <g T+" (11, 4.59), Tf (g Tf+(7, 3.96), Tf _ ®Tf +(6,3.64), 
Tf" <gTf+(5, 3.81), A 2 ~+®Tf ~ (5,3.50), Af+^Tf" (5, 3.34), Tf ~ ®Tf+(4, 3.42), 
Tf" <gTf+(4, 3.35), T x ®Tf+(4, 3.06), Tf" ® T++(4, 2.78), Tf+®A+-(4, 1.59), 
Tf-®£++(3,3.98), A+~ <g> £++(3,3.91), Tf — <g Tf+(3, 3.52), Tf" ®Tf+(3, 2.95), 
A+-®T++(3,2.56), A++ ® £+-(3,2.31), T+-&T++ (3, 2.09), Tf+®£+- (3, 1.49) 
N T++ (g)T+-(ll,4.83), T++®T+- (11, 4.63), Tf <g Tf+(8, 3.82), Tf" ® Tf+(7, 3.56), 
Tf" ®Tf+(6,3.58), A\-+®T^- (4,3.85), Af+^Tf- (4, 3.31), r+-®T++(4, 2.81), 
A+" (g £++(3,4.00), Tf" (g T 2 -+(3,3.93), Tf-®£++(3, 3.89), Tf" ®Tf +(3, 3.86), 
Tf" (g Tf +(3,3.75), T 1 —®T 2 -+ (3, 3.37), Tf" ®Tf+(3, 3.22), A+ _ ®T++(3, 3.15), 
Tf" <g T++(3,2.81), T 1 ++® J 4+-(3,2.13), A++ (g £+" (3, 2.06), Tf+®£+-(3, 1.56) 

(£,T 2 ) _ + C Tf—®T+- (12, 6.03), Tf" ®T+- (12, 5.94), Tf-^Tf" (12, 5.72), Tf" <g Tf" (12, 5.43), 
T++®Tf +(10,3.32), Tf" <g Tf" (8,4.05), Tf-^Tf" (7, 3.99), T++®Tf +(6, 3.61), 
A+-®Tf- (5,4.17), T++ ® Tf +(5,3.56), Tf" <g Tf" (4, 5.15), A^+ ® £++(4, 4.50), 
T+-®T{ "(4,4.33), A++®Tf +(4,3.93), A 2 ~+ ® £++(4, 3.76), T++ <g Tf+(4, 3.65), 
A 2 -+®T++(4,3.04), Tf" ®£+-(4,2.36), Tf" ®£+"(4, 1.88), Tf" ®A+"(4, 1.76), 
A^+®T++(4, 1.73), Tf" ®A+"(4, 1.64) 
N Tf" ®T+- (11, 5.68), Tf" ®T+- (11, 5.63), Tf " <g Tf" (11, 4.79), Tf ~ <g Tf" (10, 4.74), 
Tf" ig Tf" (7,3.61), T++®T-+ (5, 3.88), T++ <g Tf+ (5, 3.40), Tf+®Tf +(5, 2.71), 
A+-®T 2 — (4,3.91), T+-®T 2 — (4,3.33), T++ (g Tf+ (4, 2.49), ^+®T++(4, 2.17), 
Tf" ®A+"(4, 1.92), A"+ ® £++(3,4.60), T+~ (g Tf- (3, 4.42), Tf" ®£+" (3, 2.47), 
A 2 ~+ (g £++(2,4.70), T+-®T— (2,3.94), A++®Tf +(2, 3.57), A 2 -+®T++(2, 3.11), 
Tf" ®£+" (2, 2.44), Tf~ ®A+"(2, 1.63) 

A++ C Tf" ®T 2 — (7,1.25), Tf" ®T 2 ~ (7,0.89), Tf ~ (g T 2 ~ (5, 1.53), T++®T++ (2, 0.72), 

T++(gT++(l, 1.27), A+" <gA+" (1,0.50) 
N Tf-^Tf" (8,0.95), A+-®A+-{6, 1.13), Tf~ igTf" (5, 1.82), T++®T++(5, 1.15), 
Tf~ ®T 2 — (4,0.00), T++ (gT++(2, 1.77), A+" (g A+"(l, 0.30) 

Af+ C Tf" ig Tf" (6, 1.43), T+-®T 2 — (6, 0.78), T++ ig Tf +(4, 2.23), Tf+®Tf+(l, 1.09) 
N r+-®r 2 — (8, 1.18), T+- g Tf (7, 2.24), T++ g Tf +(4, 1.94), Tf+®Tf+(l, 0.97) 



ACKNOWLEDGMENTS of Canada, the Canada Research Chairs Program, the 

Canada Foundation for Innovation, and the Government 
This work was supported in part by the Natural °f Saskatchewan. 
Sciences and Engineering Research Council (NSERC) 



[1] S. Godfrey and J. Napolitano, Rev. Mod. Phys. 71, 1411 

(1999), hep-ph/9811410. 
[2] F. E. Close, Int. J. Mod. Phys. A17, 3239 (2002), hep- 

ph/0110081. 

[3] E. Klempt and A. Zaitsev, Phys. Rept. 454, 1 (2007), 

arXiv:0708.4016 [hep-ph]. 
[4] R. Frezzotti and G. C. Rossi, JHEP 08, 007 (2004), hep- 

lat/0306014. 

[5] R. Frezzotti, P. A. Grassi, S. Sint, and P. Weisz (Alpha), 

JHEP 08, 058 (2001), hep-lat/0101001. 
[6] G. M. von Hippel, R. Lewis, and R. G. Petry, Comput. 

Phys. Commun. 178, 713 (2008), arXiv:0707.2788 [hep- 

lat]. 

[7] G. M. von Hippel, R. Lewis, and R. G. Petry, PoS 
LAT2007, 043 (2007), arXiv:0710.0014 [hep-lat]. 

[8] P. Boucaud et al. (ETM) (2008), arXiv:0803.0224 [hep- 
lat]. 

[9] P. Lacock, C. Michael, P. Boyle, and P. Row- 
land (UKQCD), Phys. Rev. D54, 6997 (1996), hep- 



lat/9605025. 

[10] S. Basak et al., Phys. Rev. D72, 094506 (2005), hep- 
lat/0506029. 

[11] S. Basak et al. (Lattice Hadron Physics (LHPC)), Phys. 

Rev. D72, 074501 (2005), hep-lat/0508018. 
[12] D. Harnett, R. Lewis, and R. G. Petry, PoS LAT2006, 

194 (2006), hep-lat/0609071. 
[13] R. C. Johnson, Phys. Lett. B114, 147 (1982). 
[14] J. E. Mandula, Phys. Lett. B135, 155 (1984). 
[15] I. H. Jorysz and C. Michael, Nucl. Phys. B302, 448 

(1988). 

[16] C. W. Bernard et al. (MILC), Phys. Rev. D56, 7039 
(1997), hep-lat/9707008. 

[17] J. F. Cornwell, Group theory in physics: An introduction 
(Academic Press, San Diego, 1997). 

[18] H.-W. Lin and S. D. Cohen (2007), arXiv:0709.1902 [hep- 
lat]. 

[19] C. Michael, Nucl. Phys. B259, 58 (1985). 

[20] T. Burch, C. Gattringer, L. Y. Glozman, C. Hagen, 



23 



and C. B. Lang, Phys. Rev. D73, 017502 (2006), hep- 
lat/05 11054. 

[21] M. Luscher and U. Wolff, Nucl. Phys. B339, 222 (1990). 
[22] Available from http://www.python.org/. 
[23] K. Levenberg, Quart. Appl. Math. 2, 164 (1944). 
[24] D. Marquardt, SIAM J. Appl. math. 11, 431 (1963). 
[25] A. M. Abdel-Rehim, R. Lewis, R. M. Woloshyn, and 

J. M. S. Wu, Phys. Rev. D74, 014507 (2006), hep- 

lat/0601036. 

[26] K. Jansen, A. Shindler, C. Urbach, and I. Wetzorke 
(XLF), Phys. Lett. B586, 432 (2004), hep-lat/0312013. 

[27] A. M. Abdel-Rehim, R. Lewis, and R. M. Woloshyn, 
Phys. Rev. D71, 094505 (2005), hep-lat/0503007. 

[28] F. Farchioni et al., Nucl. Phys. Proc. Suppl. 140, 240 
(2005), hep-lat/0409098. 

[29] A. M. Abdel-Rehim, R. Lewis, R. G. Petry, and R. M. 
Woloshyn, PoS LAT2006, 164 (2006), hep-lat/0610004. 

[30] B. Efron, The Jackknife, the Bootstrap and Other Re- 
sampling Plans (SIAM, Philadelphia, 1982). 

[31] R. G. Petry, University of Regina Ph.D. thesis (2008). 

[32] W. A. Bardeen, A. Duncan, E. Eichten, N. Isgur, 



and H. Thacker, Phys. Rev. D65, 014509 (2001), hep- 
lat/0106008. 

[33] W. M. Yao et al. (Particle Data Group), J. Phys. G33, 
1 (2006). 

[34] Z.-H. Mei and X.-Q. Luo, Int. J. Mod. Phys. A18, 5713 

(2003), hep-lat/0206012. 
[35] J. N. Hedditch et al., Phys. Rev. D72, 114507 (2005), 

hep-lat/0509106. 
[36] M. S. Cook and H. R. Fiebig, Phys. Rev. D74, 094501 

(2006), hep-lat/0609010. 
[37] P. Lacock and K. Schilling (TXL), Nucl. Phys. Proc. 

Suppl. 73, 261 (1999), hep-lat/9809022. 
[38] C. Bernard et al., Phys. Rev. D68, 074505 (2003), hep- 

lat/0301024. 

[39] N. Isgur, R. Kokoski, and J. E. Paton, Phys. Rev. Lett. 
54, 869 (1985). 

[40] C. J. Morningstar and M. Peardon, Phys. Rev. D 56, 

4043 (1997). 
[41] F. E. Close (1995), hep-ph/9511442. 



